home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / cmln0286.arc / EXOTIC1.CUT < prev    next >
Text File  |  1986-02-03  |  5KB  |  112 lines
  1.  
  2.                         Sample Program #2
  3.                 From:  Exotic Language on RATFOR
  4.                           February 1986
  5.  
  6.  
  7. subroutine gaujor(a,n,detmt) #Gauss-Jordan matrix inversion with pivoting
  8. # Number of operations is n**2*(n-1) multiplies and n**2 divides
  9. # Reference:  Bauer, F. L. and Reinsch, C., "Inversion of Positive Definite
  10. #    Matrices by the Gauss-Jordan Method", given as pp 45-49 of
  11. #    Hanbook for Automatic Computation, J. H. Wilkinson and C. Reinsch, Eds.,
  12. #    Springer-Verlag, 1971.  The general case of Gauss-Jordan reduction
  13. #    is given on page 45.  Each Gauss-Jordan reduction is done by finding
  14. #    the largest element in the matrix which does not correspond to a reduced
  15. #    element.  This is the pivot element.  The operation of Gauss-Jordan
  16. #    reduction is solving for the element corresponding to the pivot column,
  17. #    then substituting for that row.  Gauss-Jordan reduction interchanges
  18. #    the rows AND columns corresponding to the pivot element.  As a result,
  19. #    both the rows and columns have to be interchanged to preserve matrix
  20. #    indexing.  The rows are switched during computation and the columns
  21. #    are switched afterward to simplify the algorithm.  The determinant is
  22. #    computed as the product of the pivot elements because the process is
  23. #    equivalent to diagonalization by Gaussian elimination (see "A First
  24. #    Course in Numerical Analysis", A. Ralston, McGraw-Hill (1965), p 400).
  25. #  Inputs:
  26. #    A(n,n)--square matrix of real elements
  27. #    n--dimension of A, limited to 30 by dimension statement
  28. #  Outputs:
  29. #    A(n,n)--inverse overstores input matrix
  30. #    detmt--determinant of input matrix
  31. dimension a(n,n),icused(30),ipivot(30,2)
  32. integer pvtrow,pvtcol
  33. detmt=1.
  34. do i=1,n                       #Reset "column used" flags
  35.   icused(i)=0
  36. do level=1,n                   #Level=no. of rows reduced
  37.   {
  38.   temp=-1.                     #Find largest element in matrix for pivot
  39.   do i=1,n
  40.     {
  41.     if(icused(i)!=0) next      #Skip columns if already used
  42.     do j=1,n
  43.       {
  44.       if(icused(j)!=0) next  #Skip rows if already used
  45.       test=abs(a(i,j))       #Use dabs in double precision versions
  46.       if(test>temp)
  47.         {
  48.         temp=test
  49.         pvtrow=i
  50.         pvtcol=j
  51.         }
  52.       }
  53.     }
  54.   pivot=a(pvtrow,pvtcol)       #Save pivot element as scalar
  55.   detmt=detmt*pivot            #Update determinant
  56.   if(temp==0.)                 #Error exit on singular matrixè    {
  57.     irank=level-1
  58.     write(3,3000) n,irank
  59.       3000 format(' ***Singular matrix in GAUJOR, order',i3,', rank',i3)
  60.     return                     #Zero detmt is error flag for calling routine
  61.     }
  62.   icused(pvtcol)=1             #Set "column used" flag
  63.   if(pvtrow!=pvtcol)           #Swap rows (if necessary) to put
  64.     {                          #  pivot element on diagonal
  65.     detmt=-detmt               #If rows are switched, toggle sign of detmt
  66.     do j=1,n
  67.       {
  68.       temp=a(pvtrow,j)
  69.       a(pvtrow,j)=a(pvtcol,j)
  70.       a(pvtcol,j)=temp
  71.       }
  72.     }
  73.   ipivot(level,1)=pvtrow       #Save indices of switched rows & cols
  74.   ipivot(level,2)=pvtcol       #  for column switching operation
  75. #Solve row pvtcol of y=A*x for x(pvtcol) and replace that row of A with
  76. #  the result.  Then, replace x(pvtcol) in the rest of the rows by
  77. #  substitution from the same result.  After doing this for all rows,
  78. #  the result is x=A^(-1)*y, and A has been replaced by A^(-1).
  79. #The signs of x(pvtcol) and y(pvtcol) are reversed to simplify the code.
  80. #The equation for row pvtcol is (using p in lieu of pvtcol for brevity):
  81. #  x(p)=[sum(j=1,n;j!=p) (a(p,j)/a(p,p))*x(j)]+(1/a(p,p))*y(p)
  82. #The equation for the other rows is:
  83. #  y(i)=[sum(j=1,n;j!=p) (a(i,j)-(a(i,p)*a(p,j)/a(p,p)))*x(i)]
  84. #       -(a(i,p)/a(p,p))*y(p)
  85. #Do the reduced row first to provide a(p,j)/a(p,p) for reducing other rows
  86.   a(pvtcol,pvtcol)=1.          #Diagonal element is special case
  87.   do j=1,n                     #Divide rest of reduced row by diagonal element
  88.     a(pvtcol,j)=a(pvtcol,j)/pivot
  89.   do i=1,n                     #Account for substituting x(pvtcol)
  90.     {                          #  into other rows
  91.     if(i==pvtcol) next         #Reduced row was done first
  92.     temp=a(i,pvtcol)           #Column for x(pvtcol) is used in whole row
  93.     a(i,pvtcol)=0.             #Accout for x(pvtcol) moved to left side
  94.     do j=1,n                   #Special case j=pvtcol OK
  95.       a(i,j)=a(i,j)-a(pvtcol,j)*temp
  96.     }
  97.   }
  98. for(level=n;level>0;level=level-1)
  99.   {                            #Inversion complete--interchange columns
  100.   j1=ipivot(level,1)           #  corresponding to pivoted rows to
  101.   j2=ipivot(level,2)           #  restore matrix indexing
  102.   if(j1==j2) next              #Skip "do nothing" swaps
  103.   do i=1,n
  104.     {
  105.     temp=a(i,j1)
  106.     a(i,j1)=a(i,j2)
  107.     a(i,j2)=temp
  108.     }
  109.   }
  110. returnèend
  111.  
  112.